function      [Sunit_opt, Skwh_opt, Scost_opt]=solve_S2( Sunit_guess,Scost_guess,Skwh_guess, ...
    params, A_t, Nbar_i, HH_bin, CollPerc_t, Pol_t, Own_t, Inc_t, ...
 P_state, P_sale_state, PIns_state, PIns0_state, nat_level, lambda, min_damages, net_cost, dW_ds_tol, ...
    Load, RenProd, PPRegion, PPRegion2, PPInt, PPCap, PPInf,PPRandStacked, ...
 ProdMedian,DamagesCoef1,DamagesCoef2, DamagesCoefPoll1,DamagesCoefPoll2, TractSunlightProfile)
    
  global nstate tract_state_xxwalk_mat  tract_nerc_xxwalk tract_state_xxwalk ntract net_metering C marginal_increase_global

  global  tl_global Sunit_tg Scost_tg Skwh_tg scc over_dist PV_hh PV trans_con_loc
    
[ sigma, gamma_0, gamma_coll, gamma_pol, gamma_own, gamma_inc, gamma_N, gamma_2N]=unpack_params(params);
    
%dW_ds_tol=.1;
dW_ds_tot=dW_ds_tol+1;

run_damages=1;

nbins=size(HH_bin,2);

%lambda

Sunit=Sunit_guess;
Scost=Scost_guess;
Skwh=Skwh_guess;

max_iter=70;

%max_iter=2
if over_dist==1
    max_iter=100;
end

if tl_global==1
    Sunit_tg= sum(tract_state_xxwalk_mat.*(ones(ntract,1)*Sunit_guess'),2); % when tract-level is on these are passes as globals
    Scost_tg= sum(tract_state_xxwalk_mat.*(ones(ntract,1)*Scost_guess'),2);
    Skwh_tg= sum(tract_state_xxwalk_mat.*(ones(ntract,1)*Skwh_guess'),2);


    Skwh_tg=Skwh_tg+Sunit_tg./A_t;
    Sunit_tg=0*Sunit_tg;
 


    max_iter=120
  %  max_iter=2
end

%Nbar_t=sum(Nbar_i.*HH_bin,2);

HH_t= sum(HH_bin,2);
HH_state=sum((HH_t*ones(1,nstate)).*tract_state_xxwalk_mat)';     
HH_nat=sum(HH_state);

HH_tc=HH_t;
HH_tc(HH_tc==0)=1;

 

P_t=sum(tract_state_xxwalk_mat.*(ones(ntract,1)*P_state'),2);
PSale_t=sum(tract_state_xxwalk_mat.*(ones(ntract,1)*P_sale_state'),2);
PIns_t=sum(tract_state_xxwalk_mat.*(ones(ntract,1)*PIns_state'),2);
PIns0_t=sum(tract_state_xxwalk_mat.*(ones(ntract,1)*PIns0_state'),2);



scale_fact_unit0=.5*ones(nstate,1);
scale_fact_cost0=.5*ones(nstate,1);
scale_fact_kwh0=2E-7*ones(nstate,1);

if scc==152 || scc==148
    scale_fact_unit0=.1*ones(nstate,1);
    scale_fact_cost0=.1*ones(nstate,1);
    scale_fact_kwh0=2E-8*ones(nstate,1);
end

if scc==330
    scale_fact_unit0=.1*ones(nstate,1);
    scale_fact_cost0=.1*ones(nstate,1);
    scale_fact_kwh0=2E-8*ones(nstate,1);
end

if tl_global==1
    scale_fact_unit0=zeros(ntract,1); % not distant from kwh subs given A constant in tract
    scale_fact_cost0=.1*ones(ntract,1);
    scale_fact_kwh0=2E-7*ones(ntract,1);




    %stopping extreme corners[
    N_mean=15*ones(ntract,1);

    % extreme corner solutions not optimal (i've checked). THis encourages
    % interior solution
    ub_sunit=1000./N_mean;
    ub_scost=1000./(PIns0_t+PIns_t.*N_mean);
    ub_skwh=1000./(N_mean.*A_t);

end

%if all hhs in state in same region, we know sunit and scost are 0
if lambda==1 && over_dist~=1
for s_i=1:nstate
   region_cut=tract_nerc_xxwalk(tract_state_xxwalk==s_i);
   
   region_min=min(region_cut);
   region_max=max(region_cut);
   
   if region_min==region_max
       scale_fact_cost0(s_i,1)=0;
       scale_fact_unit0(s_i,1)=0;
   end
    
    
end
end


if lambda==1 && over_dist~=1

    scale_fact_cost0=10*scale_fact_cost0;
end

if nat_level==1
    Sunit_nat=mean(Sunit);
    Scost_nat=mean(Scost);
    Skwh_nat=mean(Skwh);
    scale_fact_unit0=mean(scale_fact_unit0)/20;
    scale_fact_cost0=mean(scale_fact_cost0)/20;
    scale_fact_kwh0=mean(scale_fact_kwh0)/20; 
end

iter=1;

if trans_con_loc==1
    max_iter=max_iter-10;
end


while ((dW_ds_tot>dW_ds_tol) ) && iter<max_iter
           
    Sunit_old=Sunit;
    Skwh_old=Skwh;
    Scost_old=Scost;
    
    Sunit_t= sum(tract_state_xxwalk_mat.*(ones(ntract,1)*Sunit'),2);
    Scost_t= sum(tract_state_xxwalk_mat.*(ones(ntract,1)*Scost'),2);
    Skwh_t= sum(tract_state_xxwalk_mat.*(ones(ntract,1)*Skwh'),2);

    if tl_global==1

        Sunit_old_tg=Sunit_tg;
        Scost_old_tg=Scost_tg;
        Skwh_old_tg=Skwh_tg;

        Sunit_t= Sunit_tg;
        Scost_t= Scost_tg;
        Skwh_t= Skwh_tg;

       
    end
    
    N_i=solve_N_i(Nbar_i, A_t, Scost_t, Sunit_t, P_t, PSale_t,  PIns0_t, PIns_t, Skwh_t, gamma_N, gamma_2N);
    
 
    if nat_level==1
         Sunit_nat_old=Sunit_nat;
         Skwh_nat_old=Skwh_nat;
         Scost_nat_old=Scost_nat;
         Sunit=Sunit_nat*ones(nstate,1);
         Scost=Scost_nat*ones(nstate,1);
         Skwh=Skwh_nat*ones(nstate,1);
         max_iter=100;
    end
      


    [M,B, D, Dpol, MargDamagesHour,  M_bin, pi_bin, Exp_sub_t, Tot_ut_sim]=...
        simulate(params, A_t, Nbar_i, HH_bin, CollPerc_t, Pol_t, Own_t, Inc_t, ...
    Sunit, Skwh, Scost, P_state, P_sale_state, PIns_state,PIns0_state, ...
    Load, RenProd, PPRegion, PPRegion2, PPInt, PPCap, PPInf,PPRandStacked, ...
 ProdMedian,DamagesCoef1,DamagesCoef2, DamagesCoefPoll1,DamagesCoefPoll2, TractSunlightProfile, run_damages);

    pi_bin(N_i==0)=0;

    ResProd=calc_res_prod(A_t, M, TractSunlightProfile); %yearly res prod in mwh

    marg_damages_t=calc_marg_damages4(MargDamagesHour, TractSunlightProfile, ResProd); % cost in dollars of 1 unit increase in electricity in each tract
 
    dpi_ds_unit_bin=N_i.*HH_bin.*pi_bin.*(1-pi_bin)./sigma; % marginal installers wrt unit at bin level
    dN_ds_unit_bin=-(1./(2*gamma_2N));   %change in N of inframarginals  
    
    if net_metering==1       
         dN_ds_unit_bin=-(1./(2*(gamma_2N- A_t.*C.*(P_t-PSale_t))));   %change in N of inframarginals 
         dN_ds_unit_bin=dN_ds_unit_bin*ones(1,nbins);
    end   
    dN_ds_unit_bin(N_i >= Nbar_i)=0;

    dM_ds_unit_t=sum(N_i.*dpi_ds_unit_bin+dN_ds_unit_bin.*pi_bin.*HH_bin,2);
    
    dpi_ds_kwh_bin=A_t.*N_i.*HH_bin.*pi_bin.*(1-pi_bin)./sigma; % marginal installers wrt kwh at bin level



    dN_ds_kwh_bin=-(A_t./(2*gamma_2N));   %change in N of inframarginals  





    if net_metering==1
        dN_ds_kwh_bin=-(A_t./(2*(gamma_2N- A_t.*C.*(P_t-PSale_t)))); 
        dN_ds_kwh_bin=dN_ds_kwh_bin*ones(1,nbins);
    end

    if over_dist==1
        dN_ds_kwh_bin=(PV_hh./PV)*dN_ds_kwh_bin;
        dpi_ds_kwh_bin=(PV_hh./PV)*dpi_ds_kwh_bin;
    end


    dN_ds_kwh_bin(N_i >= Nbar_i)=0;
    dM_ds_kwh_t=sum(N_i.*dpi_ds_kwh_bin+dN_ds_kwh_bin.*pi_bin.*HH_bin,2); 
   
    dpi_ds_cost_bin=(PIns0_t+ PIns_t.*N_i).*HH_bin.*pi_bin.*(1-pi_bin)./sigma; % marginal installers wrt cost at bin level
    dN_ds_cost_bin=-(PIns_t./(2*gamma_2N));   %change in N of inframarginals 
    if net_metering==1
        dN_ds_cost_bin=-(PIns_t./(2*(gamma_2N- A_t.*C.*(P_t-PSale_t)))); 
        dN_ds_cost_bin=dN_ds_cost_bin*ones(1,nbins);
    end
    dN_ds_cost_bin(N_i >= Nbar_i)=0;
    dM_ds_cost_t=sum(N_i.*dpi_ds_cost_bin+dN_ds_cost_bin.*pi_bin.*HH_bin,2);    
    
    total_marg_damages_unit_t= dM_ds_unit_t.*A_t.* marg_damages_t;  % at tract level- number of additional panels (derivative) multiplied by A by marg damages    
    total_marg_damages_kwh_t= dM_ds_kwh_t.*A_t.* marg_damages_t;  % at tract level- number of additional panels  (derivative) multiplied by A by marg damages    
    total_marg_damages_cost_t= dM_ds_cost_t.*A_t.* marg_damages_t;  % at tract level- number of additional panels  (derivative) multiplied by A by marg damages    
   
    [sub_bin, sub_N_bin]=calc_sub_bin(N_i, A_t, P_t, PIns0_t, PIns_t,Sunit, Scost, Skwh);

    % government cost to marginals in each bin
    dcost_ds_unit_bin=dpi_ds_unit_bin.*sub_bin +dN_ds_unit_bin.*pi_bin .*sub_N_bin.*HH_bin;
    dcost_ds_cost_bin=dpi_ds_cost_bin.*sub_bin +dN_ds_cost_bin.*pi_bin.*sub_N_bin.*HH_bin;
    dcost_ds_kwh_bin=dpi_ds_kwh_bin.*sub_bin +dN_ds_kwh_bin.*pi_bin.*sub_N_bin.*HH_bin;
    
    %cost to marginals
    dcost_ds_unit_t=sum(dcost_ds_unit_bin,2);
    dcost_ds_cost_t=sum(dcost_ds_cost_bin,2);
    dcost_ds_kwh_t=sum(dcost_ds_kwh_bin,2);
    
    
    %cost to inframarginals
    dinfcost_ds_unit_t=sum(HH_bin.*pi_bin.*N_i,2);
    dinfcost_ds_cost_t=sum(HH_bin.*pi_bin.*(PIns0_t+ PIns_t.*N_i),2);
    dinfcost_ds_kwh_t=sum(HH_bin.*pi_bin.*A_t.*N_i,2);
    
    
    %other stuff that's useful for scaling new estimates    
    denom_unit_t=sum(N_i.*dpi_ds_unit_bin,2);
    denom_kwh_t=sum(N_i.*A_t.*dpi_ds_kwh_bin,2);   
    denom_cost_t=sum((PIns0_t+PIns_t.*N_i).*dpi_ds_cost_bin,2); 




    
    %% get to state level
    
    if nat_level==0 && tl_global==0
        
        dcost_ds_unit_state=sum((dcost_ds_unit_t*ones(1,nstate)).*tract_state_xxwalk_mat)';                 % gov cost tomarginal installers at state level
        dcost_ds_cost_state=sum((dcost_ds_cost_t*ones(1,nstate)).*tract_state_xxwalk_mat)';
        dcost_ds_kwh_state=sum((dcost_ds_kwh_t*ones(1,nstate)).*tract_state_xxwalk_mat)';

        total_marg_damages_unit_state=   sum((total_marg_damages_unit_t*ones(1,nstate)).*tract_state_xxwalk_mat)';     %damages of marginal installers at state level
        total_marg_damages_cost_state=   sum((total_marg_damages_cost_t*ones(1,nstate)).*tract_state_xxwalk_mat)';     %damages of marginal installers at state level
        total_marg_damages_kwh_state=   sum((total_marg_damages_kwh_t*ones(1,nstate)).*tract_state_xxwalk_mat)';     %damages of marginal installers at state level

        
        dinfcost_ds_unit_state=sum((dinfcost_ds_unit_t*ones(1,nstate)).*tract_state_xxwalk_mat)';                 % gov cost to inframarginal installers at state level
        dinfcost_ds_cost_state=sum((dinfcost_ds_cost_t*ones(1,nstate)).*tract_state_xxwalk_mat)';
        dinfcost_ds_kwh_state=sum((dinfcost_ds_kwh_t*ones(1,nstate)).*tract_state_xxwalk_mat)';
        

        
        dW_ds_unit=-total_marg_damages_unit_state-lambda*dcost_ds_unit_state-(lambda-1)*dinfcost_ds_unit_state;  %nstate by 1
        dW_ds_cost=-total_marg_damages_cost_state-lambda*dcost_ds_cost_state-(lambda-1)*dinfcost_ds_cost_state;  %nstate by 1
        dW_ds_kwh=-total_marg_damages_kwh_state-lambda*dcost_ds_kwh_state-(lambda-1)*dinfcost_ds_kwh_state;  %nstate by 1

        if over_dist==1
            dW_ds_kwh=-total_marg_damages_kwh_state-lambda*dcost_ds_kwh_state-(lambda-PV_hh./PV)*dinfcost_ds_kwh_state;  %nstate by 1. Because valued less by inframargianls
            
        end

        if min_damages==1
            dW_ds_unit=-total_marg_damages_unit_state-lambda*dcost_ds_unit_state-(lambda)*dinfcost_ds_unit_state;  %nstate by 1
            dW_ds_cost=-total_marg_damages_cost_state-lambda*dcost_ds_cost_state-(lambda)*dinfcost_ds_cost_state;  %nstate by 1
            dW_ds_kwh=-total_marg_damages_kwh_state-lambda*dcost_ds_kwh_state-(lambda)*dinfcost_ds_kwh_state;  %nstate by 1
        end

        if over_dist==1
            dW_ds_kwh((Skwh==0) & (dW_ds_kwh<0)) =0; % don't let per-kwh subs be negative. or else it will just set them super negative and transfer to other subs
        end


       if net_cost==1
            dW_ds_unit=-lambda*total_marg_damages_unit_state-lambda*dcost_ds_unit_state-(lambda-1)*dinfcost_ds_unit_state;  %nstate by 1
            dW_ds_cost=-lambda*total_marg_damages_cost_state-lambda*dcost_ds_cost_state-(lambda-1)*dinfcost_ds_cost_state;  %nstate by 1
            dW_ds_kwh=-lambda*total_marg_damages_kwh_state-lambda*dcost_ds_kwh_state-(lambda-1)*dinfcost_ds_kwh_state;  %nstate by 1
       end

       if marginal_increase_global==1
           marg_sub_unit=-total_marg_damages_unit_state./(dcost_ds_unit_state+dinfcost_ds_unit_state); % damages offset per marginal dollar
           marg_sub_cost=-total_marg_damages_cost_state./(dcost_ds_cost_state+dinfcost_ds_cost_state); % damages offset per marginal dollar
           marg_sub_kwh=-total_marg_damages_kwh_state./(dcost_ds_kwh_state+dinfcost_ds_kwh_state); % damages offset per marginal dollar


           mvpf_sub_unit=(-total_marg_damages_unit_state+dinfcost_ds_unit_state)./(dcost_ds_unit_state+dinfcost_ds_unit_state); % mvpf of marginal dollar
           mvpf_sub_cost=(-total_marg_damages_cost_state+dinfcost_ds_cost_state)./(dcost_ds_cost_state+dinfcost_ds_cost_state); %
           mvpf_sub_kwh=(-total_marg_damages_kwh_state+dinfcost_ds_kwh_state)./(dcost_ds_kwh_state+dinfcost_ds_kwh_state); %         

           Table=[ marg_sub_unit marg_sub_cost marg_sub_kwh  ...
               dcost_ds_unit_state dinfcost_ds_unit_state ...
               dcost_ds_cost_state dinfcost_ds_cost_state ...
               dcost_ds_kwh_state dinfcost_ds_kwh_state ...
               mvpf_sub_unit mvpf_sub_cost mvpf_sub_kwh];


           Title=[ "marg_sub_unit" "marg_sub_cost" "marg_sub_kwh"  ...
               "dcost_ds_unit_state" "dinfcost_ds_unit_state" ...
               "dcost_ds_cost_state" "dinfcost_ds_cost_state" ...
               "dcost_ds_kwh_state" "dinfcost_ds_kwh_state"...
                "mvpf_sub_unit" "mvpf_sub_cost" "mvpf_sub_kwh"];

           Table=[Title; Table]

           if over_dist==0

               title=('Results/MargSubsidy.xls');
                writematrix(Table, title); 

           end

           if over_dist==1

               title=('Results/MargSubsidyOD.xls');
                writematrix(Table, title); 

           end

           error ('Marginal Subsidies Complete') 

           stop

       end
        
        %    cost15=[ dW_ds_cost(3:5)]
        %    unit15=[ dW_ds_unit(3:5)]
        %    kwh15=[ dW_ds_kwh(3:5)]
          %  Scost15=Scost(3:5)
          %  cost15=dW_ds_cost(3:5)
       %   Skwh15=Skwh(3:5)

      %  iter
        %per dollar of revenue raised
        dW_ds_unit_tot=sum(abs(dW_ds_unit./dinfcost_ds_unit_state));
        dW_ds_cost_tot=sum(abs(dW_ds_cost./dinfcost_ds_cost_state));
        dW_ds_kwh_tot=sum(abs(dW_ds_kwh./dinfcost_ds_kwh_state));

        dW_ds_tot=dW_ds_unit_tot+dW_ds_cost_tot+dW_ds_kwh_tot
        

        
        Sunit_pass=Sunit; % these get passed out if algorithm ends (not updated ones)
        Scost_pass=Scost;
        Skwh_pass=Skwh;    
        
        %update scale fact. Only update one of the three each round
        
        
        scale_fact_unit=0;
        scale_fact_cost=0;
        scale_fact_kwh=0;
        
        if iter<=10
            scale_fact_kwh=scale_fact_kwh0;
        end

        
         if iter>10 %&& mod(iter,3)==0
            scale_fact_unit=1*scale_fact_unit0;
         end
         
         if iter>10 %&& mod(iter,3)==1
            scale_fact_kwh=1*scale_fact_kwh0;
         end
        
         if iter>10 %&& mod(iter,3)==2
            scale_fact_cost=1*scale_fact_cost0;
         end         
        
        if iter>20
            if  dW_ds_unit_tot>dW_ds_unit_tot_old
                scale_fact_unit0=scale_fact_unit0.*.9;
            end

            if  dW_ds_cost_tot>dW_ds_cost_tot_old
                scale_fact_cost0=scale_fact_cost0.*.9;
            end

            if dW_ds_kwh_tot>dW_ds_kwh_tot_old
                scale_fact_kwh0=scale_fact_kwh0.*.9;
            end    
        end
        
        dW_ds_kwh_tot_old=dW_ds_kwh_tot;
        dW_ds_cost_tot_old=dW_ds_cost_tot;
        dW_ds_unit_tot_old=dW_ds_unit_tot;


        
        Sunit=Sunit_old+dW_ds_unit.*scale_fact_unit./HH_state;
        Scost=Scost_old+dW_ds_cost.*scale_fact_cost./HH_state;
        Skwh=Skwh_old+dW_ds_kwh.*scale_fact_kwh./HH_state;        

        if over_dist==1
            Skwh(Skwh<0)=0;
        end
        
        if sum(isnan(Sunit)+isnan(Scost)+isnan(Skwh))>0
            stop
        end

        
    elseif nat_level==1 && tl_global==0

        dcost_ds_unit_nat=sum((dcost_ds_unit_t));                 % gov cost tomarginal installers at state level
        dcost_ds_cost_nat=sum((dcost_ds_cost_t));
        dcost_ds_kwh_nat=sum((dcost_ds_kwh_t));

        total_marg_damages_unit_nat=   sum((total_marg_damages_unit_t));     %damages of marginal installers at state level
        total_marg_damages_cost_nat=   sum((total_marg_damages_cost_t));     %damages of marginal installers at state level
        total_marg_damages_kwh_nat=   sum((total_marg_damages_kwh_t));     %damages of marginal installers at state level
        
        dinfcost_ds_unit_nat=sum((dinfcost_ds_unit_t));                 % gov cost to inframarginal installers at state level
        dinfcost_ds_cost_nat=sum((dinfcost_ds_cost_t));
        dinfcost_ds_kwh_nat=sum((dinfcost_ds_kwh_t));
                

        denom_unit_nat=sum((denom_unit_t));                 % gov cost tomarginal installers at state level
        denom_kwh_nat=sum((denom_kwh_t));
        denom_cost_nat=sum((denom_cost_t));
      
        %%FOC

        dW_ds_unit=-total_marg_damages_unit_nat-lambda*dcost_ds_unit_nat-(lambda-1)*dinfcost_ds_unit_nat;  
        dW_ds_cost=-total_marg_damages_cost_nat-lambda*dcost_ds_cost_nat-(lambda-1)*dinfcost_ds_cost_nat; 
        dW_ds_kwh=-total_marg_damages_kwh_nat-lambda*dcost_ds_kwh_nat-(lambda-1)*dinfcost_ds_kwh_nat;  

        if over_dist==1
            dW_ds_kwh=-total_marg_damages_kwh_nat-lambda*dcost_ds_kwh_nat-(lambda-PV_hh./PV)*dinfcost_ds_kwh_nat;  

            

        end        
        
        if min_damages==1
              dW_ds_unit=-total_marg_damages_unit_nat-lambda*dcost_ds_unit_nat-(lambda)*dinfcost_ds_unit_nat;  
              dW_ds_cost=-total_marg_damages_cost_nat-lambda*dcost_ds_cost_nat-(lambda)*dinfcost_ds_cost_nat; 
              dW_ds_kwh=-total_marg_damages_kwh_nat-lambda*dcost_ds_kwh_nat-(lambda)*dinfcost_ds_kwh_nat;  
        end
        
        if net_cost==1
          dW_ds_unit=-lambda*total_marg_damages_unit_nat-lambda*dcost_ds_unit_nat-(lambda-1)*dinfcost_ds_unit_nat;  
          dW_ds_cost=-lambda*total_marg_damages_cost_nat-lambda*dcost_ds_cost_nat-(lambda-1)*dinfcost_ds_cost_nat; 
          dW_ds_kwh=-lambda*total_marg_damages_kwh_nat-lambda*dcost_ds_kwh_nat-(lambda-1)*dinfcost_ds_kwh_nat;  
       
        end
        

        %iter;
        %per unit of revenue raise
        dW_ds_unit_tot=sum(abs(dW_ds_unit./dinfcost_ds_unit_nat));
        dW_ds_cost_tot=sum(abs(dW_ds_cost./dinfcost_ds_cost_nat));
        dW_ds_kwh_tot=sum(abs(dW_ds_kwh./dinfcost_ds_kwh_nat));


        dW_ds_tot=dW_ds_unit_tot+dW_ds_cost_tot+dW_ds_kwh_tot

        Sunit_pass=Sunit; % these get passed out if algorithm ends (not updated ones)
        Scost_pass=Scost;
        Skwh_pass=Skwh;    
        
        %update scale fact. Only update one of the three each round
        
        
        scale_fact_unit=0;
        scale_fact_cost=0;
        scale_fact_kwh=0;
        
        if iter<10
            scale_fact_kwh=scale_fact_kwh0;
        end

        
         if iter>=10 %&& mod(iter,3)==0
            scale_fact_unit=1*scale_fact_unit0;
         end
         
         if iter>=10 %&& mod(iter,3)==1
            scale_fact_kwh=1*scale_fact_kwh0;
         end
        
         if iter>=10 % && mod(iter,3)==2
            scale_fact_cost=1*scale_fact_cost0;
         end 
         
         if iter>50 && mod(iter,3)==0
            scale_fact_unit=.1*scale_fact_unit0;
         end
         
         if iter>50 && mod(iter,3)==1
            scale_fact_kwh=.1*scale_fact_kwh0;
         end
        
         if iter>50 && mod(iter,3)==2
            scale_fact_cost=.1*scale_fact_cost0;
         end          
        

 
        Sunit_nat=Sunit_nat_old+dW_ds_unit.*scale_fact_unit./HH_nat;
        Scost_nat=Scost_nat_old+dW_ds_cost.*scale_fact_cost./HH_nat;
        Skwh_nat=Skwh_nat_old+dW_ds_kwh.*scale_fact_kwh./HH_nat;
         
        
        %cannot serve as tax
       % S(S<0)=0;
        
        if sum(isnan(Sunit_nat)+isnan(Scost_nat)+isnan(Skwh_nat))>0
            stop
        end        
    
        
    end 

    if  tl_global==1

        
        dW_ds_unit=-total_marg_damages_unit_t-lambda*dcost_ds_unit_t-(lambda-1)*dinfcost_ds_unit_t;  %nstate by 1
        dW_ds_cost=-total_marg_damages_cost_t-lambda*dcost_ds_cost_t-(lambda-1)*dinfcost_ds_cost_t;  %nstate by 1
        dW_ds_kwh=-total_marg_damages_kwh_t-lambda*dcost_ds_kwh_t-(lambda-1)*dinfcost_ds_kwh_t;  %nstate by 1

        if over_dist==1

            dW_ds_kwh=-total_marg_damages_kwh_t-lambda*dcost_ds_kwh_t-(lambda-PV_hh./PV)*dinfcost_ds_kwh_t;  
        end   


        if min_damages==1
            dW_ds_unit=-total_marg_damages_unit_t-lambda*dcost_ds_unit_t-(lambda)*dinfcost_ds_unit_t;  %nstate by 1
            dW_ds_cost=-total_marg_damages_cost_t-lambda*dcost_ds_cost_t-(lambda)*dinfcost_ds_cost_t;  %nstate by 1
            dW_ds_kwh=-total_marg_damages_kwh_t-lambda*dcost_ds_kwh_t-(lambda)*dinfcost_ds_kwh_t;  %nstate by 1
        end

        dW_ds_unit=0*dW_ds_unit;

        pi_mean=mean(pi_bin,2);

        dW_ds_unit(pi_mean>.9)=0;
        dW_ds_cost(pi_mean>.9)=0;
        dW_ds_kwh(pi_mean>.9)=0;

        dW_ds_unit(pi_mean==0)=0;
        dW_ds_cost(pi_mean==0)=0;
        dW_ds_kwh(pi_mean==0)=0;


%non negativity constraint
        dW_ds_unit((Sunit_tg==0) & (dW_ds_unit<0)) =0;
        dW_ds_cost((Scost_tg==0) & (dW_ds_cost<0)) =0;
        dW_ds_kwh((Skwh_tg==0) & (dW_ds_kwh<0)) =0;

     

        dW_ds_unit((Sunit_tg>ub_sunit) & (dW_ds_unit>0)) =0;
        dW_ds_cost((Scost_tg>ub_scost) & (dW_ds_cost>0)) =0;
        dW_ds_kwh((Skwh_tg>ub_skwh) & (dW_ds_kwh>0)) =0;



        % some places have no possible installers so dividing by 0 here.
        dinfcost_ds_unit_t(dinfcost_ds_unit_t==0)=1;
        dinfcost_ds_cost_t(dinfcost_ds_cost_t==0)=1;
        dinfcost_ds_kwh_t(dinfcost_ds_kwh_t==0)=1;
    
        dW_ds_unit_tot=sum(abs(dW_ds_unit./dinfcost_ds_unit_t));
        dW_ds_cost_tot=sum(abs(dW_ds_cost./dinfcost_ds_cost_t));
        dW_ds_kwh_tot=sum(abs(dW_ds_kwh./dinfcost_ds_kwh_t));



        dW_ds_tot=dW_ds_unit_tot+dW_ds_cost_tot+dW_ds_kwh_tot;
        
      %  kwh_t_diff=dW_ds_kwh(1:3)'
     %   cost_t_diff=dW_ds_cost(1:3)'
        
        %update scale fact. Only update one of the three each round
        
        
        scale_fact_unit=0;
        scale_fact_cost=0;
        scale_fact_kwh=0;



      %  scale_fact_unit0(Sunit_tg>ub_sunit)=scale_fact_unit0(Sunit_tg>ub_sunit)./5;
      if iter>1
        scale_fact_cost0((abs(dW_ds_cost)>abs(dW_ds_cost_old)) & (sign(dW_ds_cost)~=sign(dW_ds_cost_old)))=...
            scale_fact_cost0((abs(dW_ds_cost)>abs(dW_ds_cost_old)) & (sign(dW_ds_cost)~=sign(dW_ds_cost_old)))./2;
        scale_fact_kwh0((abs(dW_ds_kwh)>abs(dW_ds_kwh_old)) & (sign(dW_ds_kwh)~=sign(dW_ds_kwh_old)))...
            =scale_fact_kwh0((abs(dW_ds_kwh)>abs(dW_ds_kwh_old)) & (sign(dW_ds_kwh)~=sign(dW_ds_kwh_old)))./2;
      end

        dW_ds_cost_old=dW_ds_cost;
        dW_ds_kwh_old=dW_ds_kwh;
        
        if iter<=10
            scale_fact_kwh=scale_fact_kwh0;
        end

        

         
         if iter>10 %&& mod(iter,3)==1
            scale_fact_kwh=1*scale_fact_kwh0;
         end
        
         if iter>20 %&& mod(iter,3)==2
            scale_fact_cost=1*scale_fact_cost0;
         end         
        
        if iter>20

            if  dW_ds_cost_tot>dW_ds_cost_tot_old
                scale_fact_cost0=scale_fact_cost0.*.9;
            end

            if dW_ds_kwh_tot>dW_ds_kwh_tot_old
                scale_fact_kwh0=scale_fact_kwh0.*.9;
            end    
        end
        
        dW_ds_kwh_tot_old=dW_ds_kwh_tot;
        dW_ds_cost_tot_old=dW_ds_cost_tot;
        dW_ds_unit_tot_old=dW_ds_unit_tot;




        Sunit_pass=Sunit_tg; 
        Scost_pass=Scost_tg;
        Skwh_pass=Skwh_tg;  

       
       % Sunit_tg=0*Sunit_old_tg+dW_ds_unit.*scale_fact_unit./HH_tc; % keep
       % at 0
        Scost_tg=Scost_old_tg+dW_ds_cost.*scale_fact_cost./HH_tc;
        Skwh_tg=Skwh_old_tg+dW_ds_kwh.*scale_fact_kwh./HH_tc;     

       


        %can't serve as tax
        scale_fact_unit0(Sunit_tg<0)=scale_fact_unit0(Sunit_tg<0)./5;
        scale_fact_cost0(Scost_tg<0)=scale_fact_cost0(Scost_tg<0)./5;
        scale_fact_kwh0(Skwh_tg<0)=scale_fact_kwh0(Skwh_tg<0)./5;


        Sunit_tg(Sunit_tg<0)=0;
        Scost_tg(Scost_tg<0)=0;
        Skwh_tg(Skwh_tg<0)=0;


        %upper bound
        scale_fact_unit0(Sunit_tg>ub_sunit)=scale_fact_unit0(Sunit_tg>ub_sunit)./5;
        scale_fact_cost0(Scost_tg>ub_scost)=scale_fact_cost0(Scost_tg>ub_scost)./5;
        scale_fact_kwh0(Skwh_tg>ub_skwh)=scale_fact_kwh0(Skwh_tg>ub_skwh)./5;

        Sunit_tg(Sunit_tg>ub_sunit)=ub_sunit(Sunit_tg>ub_sunit);
        Scost_tg(Scost_tg>ub_scost)=ub_scost(Scost_tg>ub_scost);
        Skwh_tg(Skwh_tg>ub_skwh)=ub_skwh(Skwh_tg>ub_skwh);

        
        %cannot serve as tax
       % S(S<0)=0;
        
        if sum(isnan(Sunit_tg)+isnan(Scost_tg)+isnan(Skwh_tg))>0
            stop
        end


    end
    
    iter=iter+1;
end

Sunit_opt=Sunit_pass;
Scost_opt=Scost_pass;
Skwh_opt=Skwh_pass;   


if tl_global==1

        Sunit_tg=Sunit_pass; 
        Scost_tg=Scost_pass;
        Skwh_tg=Skwh_pass;  

            % these are just so code can run. Average subsidy at state level. Not used for simulations just for results tables
        Sunit_opt=sum(((HH_t.*Sunit_tg)*ones(1,nstate)).*tract_state_xxwalk_mat)'./(HH_state); 
        Scost_opt=sum(((HH_t.*Scost_tg)*ones(1,nstate)).*tract_state_xxwalk_mat)'./(HH_state); 
        Skwh_opt=sum(((HH_t.*Skwh_tg)*ones(1,nstate)).*tract_state_xxwalk_mat)'./(HH_state); 
end


dW_ds_tot